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A Cl VIRTUAL ELEMENT METHOD FOR THE CAHN-HILLIARD 
EQUATION WITH POLYGONAL MESHES 

PAOLA F. ANTONIETTI*, LOURENCO BEIRAO DA VEIGAt, SIMONE SCACCHI t, 

AND MARCO VERANI § 

Abstract. In this paper we develop an evolution of the virtual elements of minimal degree 
for the approximation of the Cahn-Hilliard equation. The proposed method has the advantage of 
being conforming in and making use of a very simple set of degrees of freedom, namely 3 degrees 
of freedom per vertex of the mesh. Moreover, although the present method is new also on triangles, 
it can make use of general polygonal meshes. As a theoretical and practical support, we prove the 
convergence of the semi-discrete scheme and investigate the performance of the fully discrete scheme 
through a set of numerical tests. 
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1. Introduction. The study of the evolution of transition interfaces, which is 
of paramount importance in many physical/biological phenomena and industrial pro¬ 
cesses, can be grouped into two macro classes, each one corresponding to a different 
method of dealing with the moving free-boundary: the sharp interface method and 
the phase-field method. In the sharp interface approach, the free boundary is to be 
determined together with the solution of suitable partial differential equations where 
proper jump relations have to imposed across the free boundary. In the phase field 
approach, the interface is specified as the level set of a smooth continuos function 
exhibiting large gradients across the interface. 

Phase field models, which date back to the works of Korteweg m, Cahn and 
Hilliard [T3j [301 1311 ? Landau and Ginzburg [3l] and van der Waals [33] , have been 
classicaly employed to describe phase separation in binary alloys. However, recently 
Cahn-Hilliard type equations have been extensively used in an impressive variety of 
applied problems, such as, among the others, tumor growth imisn], origin of Saturn’s 
rings iia, separation of di-block copolymers m, population dynamics HZ], image 
processing [9] and even clustering of mussels m- 

Due to the wide spectrum of applications, the study of efficient numerical methods 
for the approximate solution of the Cahn-Hilliard equation has been the object of an 
intensive research activity. Summarizing the achievements in this field is a tremendous 
task that go beyond the scope of this paper. Here, we limit ourvselves to some 
remarks on finite element based methods, as the main properties (and limitations) of 
these schemes are instrumental to motivate the introduction of our new approach. As 
the Cahn-Hilliard equation is a fourth order nonlinear problem, a natural approach 
is the use of finite elements (FEM) as in [23 [21]. However, in order to avoid 
the well known difficulty met in the implementation of finite elements, another 
possibility is the use of non-conforming (see, e.g., [22]) or discontinuous (see, e.g., 
ESI) methods; the drawback is that in such case the discrete solution will not satisfy 
a regularity. Alternatively, the most common strategy employed in practice to 
solve the Cahn-Hilliard equation with (continuos and discontinuous) finite elements is 
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to use mixed methods (see e.g. [23j [23] and [32] for the continuous and discontinuous 
setting, respectively). Clearly, the drawback of this approach is the increase of the 
numbers of degrees of freedom, and thus of the computational cost. Very recently, the 
difficulty related to the practical use of basis functions has been addressed with 
success also in the framework of isogeometric analysis m- 

In this paper, we introduce and analyze the virtual element method (VEM) 
for the approximate solution of the Cahn-Hilliard equation. This newly introduced 
method (see, e.g., [ 3 ] for an introduction to the method and [ 6 ] for the details of 
its practical implementation) is characterized by the capability of dealing with very 
general polygonal/polyedral meshes and to possibility of easily implementing highly 
regular discrete spaces. Indeed, by avoiding the explicit construction of the local 
basis functions, the VEM can easily handle general polygons/polyhedrons without 
complex integrations on the element. In addition, thanks to this added flexibility, it 
was discovered [HE] that virtual elements can also be used to build global discrete 
spaces of arbitrary regularity {C^ and more) that are quite simple in terms of degrees 
of freedom and coding. Other virtual element contributions are, for instance Una 
ElIHlIia EZl El EZj, while for a very short sample of other EEM-inspired methods 
dealing with general polygons we refer to [TOj [161 HI UHl EEl SHI SH SS] • 

In the present contribution we develop a modification of the virtual elements 
(of minimal degree) of [7] for the approximation of the Cahn-Hilliard equation. Also 
taking inspiration from the enhancement techniques of [ 2 ] , we define the virtual space 
in order to be able to compute three different projection operators, that are used 
for the construction of the discrete scheme. Afterwards, we prove the convergence of 
the semi-discrete scheme and investigate the performance of the fully discrete scheme 
numerically. We underline that, on our knowledge, this is the first application of the 
newborn virtual element technology to a nonlinear problem. 


The paper is organized as follows. In Section]^ we describe the proposed virtual 
element method. In Section]^ we develop the theoretical error estimates. In Section]^ 
we present the numerical tests. 


2. The continuous and discrete problems. In this section, after presenting 
the Cahn-Hilliard equation, we introduce the Virtual Element discretization. The 
proposed strategy takes the steps from the methods described in ram for the 
Kirchhoff and Poisson problems, respectively, combined with an enhancement strat¬ 
egy first introduced in [2]. The present virtual scheme makes use of three different 
projectors and of a particular construction to take care of the nonlinear part of the 
problem. 

2.1. The continuous problem. Let C be an open bounded domain. Let 
2 p{x) = (1 — x^)^/4 and let 0(x) = '^'(x), we consider the following Cahn-Hilliard 
problem: find u{x,t) : O x [0,T] ^ M such that: 


( 2 . 1 ) 


' dfU — A(0(14) — = 0 

< 'u(-,0) = 'Uo(-) 

^ dn'u = 5n(0('u) - j‘^Au{t)) = 0 


in X [0, T], 
in f), 

on dfl X [0, T], 


where 5n denotes the (outward) normal derivative and 7 G 0 < 7 <C 1, represents 
the interface parameter. Throughout the paper we will employ the standard notation 
for Sobolev spaces [T]. We now introduce the variational form of (2.1) that will be 
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used to derive the virtual element discretization. To this aim, we preliminary define 
the following bilinear forms 

a^{v^w) = f (V^u) : (V^re) dx Vu, le G 

Jn 

a^{v,w)= f (Vv) ' {\/w) dx Vu, le G 

Jn 

a°{v,w) = [ 

Jn 


vwdx 


Wv,w G L\n), 


and the semi-linear form 


r{z;v,w) = / (j)'{z)\/V ' w dx ^z^v^w ^ H^{Q) 

Jn 

where all the symbols above follow a standard notation. Finally, introducing the space 

(2.2) V = [v ^ : dnU = 0 on 

the weak formulation of problem reads as: find G V such that 


(2.3) 


oP{dtU, v) + 7^a^(i4, v) + r{u] u,v) = d Vu G V, 
'u(-,0) = 'Uo(-)- 


In the theoretical analysis of Section we will work under the following regularity 
assumption on the solution of ( |2.3[ ) 

(2.4) 

see, e.g., m for a possible proof under higher regularity hypotheses on the initial 
datum uq. 

2.2. A Virtual Element space. In the present section we describe the 
virtual element space Wh C that we will use in the next section to build a 

discretization of problem ( |2.3| ). From now on, we will assume that ^2 is a polygonal 
domain in M^. 

Our construction will need a few steps. Let Clh represent a decomposition of ^2 
into general polygonal elements E of diameter He- In the following, we will denote 
by e the edges of the mesh and, for all e G dE^ will denote the unit normal 
vector to e pointing outward to E. We will use the symbol P/c(cl;) to denote the space 
of polynomials of degree less than or equal to k living on the set cj C Finally, we 
will make use of the following local bilinear forms for all E £ 

a^(v,w)= [ : (V^w) drc 'iv,w e 

Je 

a^{v,w) = f (Vu) • (Vie) dx Vu, le G iL^(F^), 

Je 

a%{v,w) = / vwdx \/v,w e L‘^{E). 

Je 


(2.5) 
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Given an element E e the augmented local space Vh\E is defined by 
( 2 ^ 6 ) 

Vh\E = {t' e H\E) ; e P2(E), v\aE e C\dE),v\, G P3(e) Ve G dE, 

Vv\aE e [C\dE)f,d^v\, G Pi(e) Ve G Se], 


with 5n denoting the (outward) normal derivative. The space Vh\E is made of functions 
that are continuous and piecewise cubic on the boundary, with continuous gradient 
on the boundary, normal linear component on each edge and such that its bilaplacian 
is a quadratic polynomial. _ 

We now introduce two sets D1 and D2 of linear operators from Vh\E into M. For 
all Vh G Vh\E they are defined as follows: 

D1 contains linear operators evaluating at the n = n{E) vertexes of E; 

D2 contains linear operators evaluating \/vh at the n = n{E) vertexes of E. 
Note that, as a consequence of definition (2.6), the output values of the two sets of 
operators D1 and D2 are sufficient to uniquely determine Vh and \/vh on the boundary 
of E. 

Let us now introduce the projection operator 11^: Vh\E ^ 2(^)5 defined by 


i^E^h,q) =aE{vh,q) Vg G P2 (E) 
{{n^Vh,q))E = {{vh,q))E VgGPi(E), 


for all Vh G Vh\E where ((’,-))b represents an euclidean scalar product acting on the 
function vertex values, i.e. 


{{Vh,Wh))E= y] Vh{v)Wh{v) \/Vh,Wh e C°{E). 

u vertexes 

of dE 


Some explanation is in order to motivate the construction of the operator 11^. First, 
we note that the bilinear form a^(-, •) has a non-trivial kernel, given by Pi(F^). Hence, 
the role of the second condition in ( |2.7[ ) is to select an element of the kernel of the 
operator. Moreover, it is easy to check that the operator H^ is well defined, as for 
all Vh G Vh\E h returns one (and only one) function G IP 2 (^)- Second, it is 

crucial to remark that the operator H^ is uniquely determined on the basis of the 
informations carried by the linear operators in D1 and D2. Indeed, it is sufficient to 
perform a double integration by parts on the right hand side of ( |2.7[ ), which gives 

aE{vh,q) = [ : V^qdx = [ (V^(g)n£;) ■ Vi;?,ds - [ Vh{divV‘^q) ■ 11%, 

JE JdE JdE 


and to observe that the above term on the right hand side only depends on the 
boundary values of Vh and Vvh- We note that the same holds for the right hand side 
of the second equation in (2.7), since it depends only on the vertex values of Vh- To 


conclude, as for any Vh G Vh\E^ Ihe output values of the linear operators in Dl and 
D2 are sufficient to define Vh and \/vh on the boundary, it turns out that the operator 
is uniquely determined on the basis of the evaluations performed by the linear 
operators in and D2. 

We are now ready to define our virtual local spaces 


( 2 . 8 ) Wh\E — 


h\E 




L 


'nEivh)qdx= I Vhqdx Vg G P2(-B)}- 
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We observe that, since Wh\E C Vh\E^ the operator 11^ is well defined on Wh\E and 
computable only on the basis of the output values of the operators in D1 and D2. 

Moreover, we have the following result. 

Lemma 2.1. The set of operators D1 and D2 eonstitutes a set of degrees of 
freedom for the spaee Wh\E‘ 

Proof We start by noting that the space Vh\E is associated to a well posed 
biharmonic problem on E with Dirichlet boundary data and standard volume loading, 
i.e.. 


— A^Vh assigned in E, 

Vh and dnVh assigned on dE. 


Thus the dimension of V^^e equals the dimension of the data space (loading and 
boundary data spaces). We now recall that, as already noted, the operators D1 and 
D2 uniquely determine Vh and \/vh on the boundary of E and thus the cardinality 
^{Dl} + ^{D2} exactly corresponds to the dimension of the boundary data in the 
above biharmonic problem. Therefore, since the loading data space has dimension 
equal to dim(P 2 (^)), we have 

dim(y;,|^) = #{1^1} + #{D2} + dim(P2(^)). 


Now, we observe that the space Wh\E is a subspace of Vh\E obtained by enforcing 
the constraints in (2.8), i.e a set of n linear equations, with n = dim(P 2 (^))- Since 
such equations could, in principle, not be linearly independent, all we can say on the 
dimension of Wh\E is 


(2.9) 


dimiW^E) > dim{V^E) - dim(P2(^)) = #{D1} + #{D2}. 


The proof is therefore complete if we show that any G Wh\E that vanishes on 
D1 and D2 is indeed the zero element of Wh\E- Let Vh G Wh\E vanish on HI and 
D2. First of all, this easily implies that Vh and \/vh are null on the boundary dE. 
Moreover, since the operator 11^ is linear and depends only on the output values of 


the operators in D1 and H2, it must hold = 0. Recalling definition (2.8), this 

in turn yields 


( 2 . 10 ) 


/ Vhqdx = 0 

Je 


V^GP2(^). 


Since Vh G Wh\E ^ Vh\E i we have A‘^Vh G P 2 (^)- Therefore, we can take q = A‘^Vh as 
a test function in (2.10). A double integration by parts, using also that and \/vh 


are null on 9H, then gives 


/ Vh A^Vh dx= Avh Avh dx. 
Je Je 


Thus Avh = 0 and the proof is complete by recalling again the boundary conditions 
on Vh- □ 

The space Wh\E satisfies also the following properties. The first one is that 

^2{E) C Wh\E. 
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that will guarantee the good approximation proj^rties for the space. The above 
inclusion is easy to verify, since clearly P 2 (^) ^ and the additional condition 

in ( |2.8[ ) is satisfied by P 2 (^) polynomials (being 11^ a projection on such polynomial 
space). The second property is that also the standard projection operator 11^ : 
Wh\E ^ 2 {E) is computable (only) on the basis of the values of the degrees of 
freedom D1 and D2. Indeed, for all fhe function Il%Vh G IP 2 (^) is defined 

by 

(2.11) a%{U%Vh,q) = a%{vh,q) Vg G P2(£'), 

where the bilinear form , introduced in ( |2.5[ ), represents the E{E) scalar 

product. Due to the particular property appearing in definition (2.8), the right hand 
side in (2.11) is computable using and thus U.%Vh depends only on the values of 

the degrees of freedom Dl, D2 attained by Vh and Actually, it is easy to check 

that on the space Wh\E the projectors 11^ and 11^ are the same operator (although 
for the sake of clarity we prefer to keep the notation different). 

We introduce an additional projection operator that we will need in the following. 
Wh\E ^ 2 {E) by 


We define n| 


( 2 . 12 ) 


( a^{U^Vh,q) = a^{vh,q) G P2(^) 

I / U^Vh dx = Vh dx. 

V J E J E 


We remark that, since the bilinear form a^(-, •) has a non trivial kernel (given by the 
constant functions) we added a second condition in order to keep the operator 11^ 


well defined. It is easy to check that the right hand side in (2.12) is computable on 
the basis of the values of the degrees of freedom D1 and D2. For the first equation 
in (2.12), this can be shown with an integration by parts (similarly as already done 
for the projector) 

/ Vvh'Vqdx =-{/^q)\E / Vhdx ^ / Vhdr^qds 
JE JE JdE 


and noting that the identity 


/ Vhdx = / IV%Vhdx 
J E J E 


allows to compute the integral of Vh on E using only the values of the degrees of free¬ 
dom D1 and D2. For the ease of the reader, we summarize what we have accomplished 
so far in the following remark. 

Remark 2.1. We have introdueed a set of loeal spaees Wh\E (well defined on 
general polygons and eontaining P 2 {E)) and the assoeiated loeal degrees of freedom. 
We have moreover shown that we have three different projeetion operators (eaeh one 
assoeiated to a different bilinear form appearing in the problem) that ean be eomputed 
making use only of the values of sueh degrees of freedom. 

The global discrete space can now be assembled in the classical finite element 
fashion, yielding 

Wh = {veV : v\Ee WhiE G ^^4- 


Note that, by gluing in the standard way the degrees of freedom, the ensuing functions 
will have continuous values and continuous gradients across edges. Therefore the 
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resulting space is indeed contained in and will yield a conforming solution. 

The global degrees of freedom will simply be 

• Evaluation of Vh at the vertexes of the mesh 

• Evaluation of \/vh at the vertexes of the mesh 

Thus the dimension of Wh is three times the number of vertexes in the mesh. As 
a final note we observe that, in practice, it is recommended to scale the degrees of 
freedom D2 by some local characteristic mesh size in order to obtain a better 
condition number of the final system. 


2.3. Virtual forms. The second key step in the contruction of the method is the 
definition of suitable discrete forms. Analogously to the finite element case, this forms 
will be constructed element by element and will depend on the degrees of freedom of 
the discrete space. Unlike in the finite element case, this forms will not be obtained by 
some Gauss integration of the shape functions (that are unknown inside the elements) 
but rather using the projection operators that we defined in the previous section. 

We start by introducing a discrete approximation of the three exact local forms 
in (2.5). By making use of the projection operators of the previous section, the 
development of the bilinear forms follows a standard approach in the virtual element 
literature. We therefore refer, for instance, to [3] for more details and motivations 
regarding this construction. Let E G be any element of the polygonal partition. 
We introduce the following (strictly) positive definite bilinear form on Wh\E x Wh\E 


SEiVh,Wh) = y] (vhiu) Wh{u) + (Kf Vvhiu) ■Vwh{u)^ yvh,Wh€Wh\E, 

ly vertexes 

of dE 


where hj^ is some characteristic mesh size lenght associated to the node u (for instance 
the maximum diameter among the elements having as a vertex). 

Recalling (2.5), we then propose the following discrete (and symmetric) local 
forms 


+ h^‘^SE{Vh ~ ~ n^UZ/j,), 

(2.13) aZsi^h, Wh) = a%{n'^Vh, U%Wh) + SE{vh - W^Vh, Wh - U%Wh), 

'^h) = U^Wh) + SE(Vh ~ ^%Vh, Wh — Il%Wh), 


for all Vh,WhG Wh\E- 

The consistency of the discrete bilinear forms is assured by the first term on the 
right hand side of each relation, while the role of the second term se(’, ’) is only to 
guarantee the correct coercivity properties. Indeed, noting that the projection oper¬ 
ators appearing above are always orthogonal with respect to the associated bilinear 
form, it is immediate to check the following consistency lemma. 

Lemma 2.2 (consistency). For all the three bilinear forms in (2.13) it holds 


^h) = at (fo '>^h) Vp e F2{E), \/Vh e Wh\E, 


where the symbol f stands for the symbol A, V or 0. 

The lemma above states that the bilinear forms are exact whenever one of the 
two entries is a polynomial in P 2 (^)- In order to present a stability result for the 
proposed discrete bilinear forms, we need some mesh regularity assumptions on the 
mesh sequence 





Assumption 2.1. We assume that there exist positive eonstants Cg and c'^ sueh 
that every element E G is star shaped with respeet to a ball with radius p > CghE 

and every edge e G dE has at least length he > c'^hE- 

Under the above mesh regularity conditions, we can show the following lemma. 
Since the proof is standard and based on a scaling argument, it is omitted. 

Lemma 2.3 (stability). Let Assumption 2.1 hold. There exist two positive eon¬ 


stants c*,c* independent of the element E G sueh that 

c*a^Ei^h,Vh) < al E{vh,Vh) < I'/i) G Wh\E, 

where the symbol f stands for the symbol A, V or 0. 

Note that, as a consequence of the above lemma, it is immediate to check that 
the bilinear forms £;(*, •) continuous with respect to the relevant norm: for 


(2.13)^, for ( 2 . 13)2 ( 2 . 13 ) 3 . global discrete bilinear forms will be 

written (following the classical finite element procedure) 


ai{vh,Wh)= al^^{vh,Wh) Vvh,WheWh, 

with the usual multiple meaning of the symbol f. 

We now turn our attention to the semilinear form r(-;-,-), that we here write 
more explicitly: 


r{z;v^w) = rE{z;v,w) Mz^v^w ^ 

Eefih 

rE{z;v,w) = / {3z{x)‘^ — l)\/v{x) • \/w{x) dx MEe^h- 
JE 

On each element E^ we approximate the term w{x)^ with its average, computed using 
the L^{E) bilinear form £;(*, •)• 

{wI)\e ^ \E\-^al E{wh,Wh), 

where \E\ denotes the area of element E. This approach will turn out to have the 
correct approximation properties and, moreover, it preserves the positivity of w‘^. We 
therefore propose the following approximation of the local nonlinear forms 

rh,EiZh;Vh,Wh) = (l)'{Zh)\E 0'h,Ei^h,Wh) '^Zh,Vh,Wh € Wh\E 

where (l)'{zh)^E~^\^\~^^h,Ei^h^ Zh) — 1. The global form is then assembled as usual 


rh{zh;vh,Wh) = rh,E{zh\Vh,Wh) 'iwh,rh,Vh eWh- 

Eedh 


2.4. Discrete problem. We here outline the Virtual Element discretization of 
problem (2.3), that follows a Galerkin approach in space combined with a backward 
Euler in time. Let us introduce the space with boundary conditions 


= WhnV = {v eWh : = 0 on dQ}. 
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As usual, it is convenient to first introduce the semi-discrete problem: 

(2.14) 

{ Find Uh{',t) in such that 

al{dtUh, Vh) + j‘^a^{uh, Vh) + rh{uh, uh] vh) = 0 Vu/, G VF^, a.e. in (0, T), 

with ^ a suitable approximation of and where the discrete forms above 
have been introduced in the previous section. 


In order to introduce the fully discrete problem, we subdivide the time interval 
[0, T] into N uniform sub-intervals of length k = T/N by selecting, as usual, the time 
nodes 0 = to < < ... < tjsf-i <t]sf=T. We now search for with 

G representing the solution at time 

The fully discrete problem reads as follows: Given ^ for ^ = 

1,..., At look for G bF^ such that 

(2.15) k-'^al{vi^ - ,Vh) + '^'^a^{vi^,Vh) + rh{vik,vik\Vh) =0 Mvu & W^. 


3. Error analysis of the semi-discretization scheme. Throughout the sub¬ 
sequent discussion, we will employ the notation x <yto denote the inequality x < Cy 
being C a positive constant independent of the discretization parameters but that may 
depend on the regularity of the underlying continuous solution. Moreover, note that 
(unless needed to avoid confusion) in the sequel the dependence of u and Uh on time 
t is left implicit and the bounds involving u or hold for all t G (0,T]. 

In this section we present the convergence analysis of the semidiscrete Virtual 


Element formulation given in (2.14). Our theoretical analysis will deal only with the 


semi-discrete case since the main novelty of the present paper is the (virtual element) 
space discretization. The error analysis of the fully discrete scheme follows from the 
analysis of the semi-discrete case employing standard techniques as for in the classical 
finite element case (see, e.g, [iT]b 

The subsequent convergence analysis will be performed under the following well 
accepted regularity assumption on the semi-discrete solution of (2.14) (see, e.g., 
[22] for a discussion on its validity). 

Assumption 3.1. The solution Uh of (2.14) satisfies 


Uh G L^{n) Vt G (0,T]. 


As a starting point, we recall the following approximation 


[371111. 

Proposition 3.1. 
H^{E) there exists G 


Assume that Assumption 2A is satised. 
¥k{E), k > 0 and vj G Wh\E ^ueh that 


result, see m and 
Then for every v G 


1'^ ^ \'^\ h ^{ e ) ’ 1 < s < A: H- 1, — 0,1,..., 5, 

1^ ~ ^ s = 2,3, ^ = 0,1,..., s. 


where the hidden eonstant depends only on k and on the eonstants in Assumption \2. 1\ 
Let 


(j)'{u)^E=^X\ ^a%{u,u)-l 
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we define 


rh(u-,Vh,Wh)= ^ (j>'{u)^j^ajE(yh,Wh). 

We introduce the elliptic projection P^v G for v G defined by 

(3.2) bh{P^v, xph) = ( 7 ^ - V ■ {(p'{u)Vv) + av, xph) 

for all ^ where bh{'r) is the bilinear form 

(3.3) bh{vh, Wh) = 7 ^a^(t>ft, wn) + fh{u-, vu, wu) + a(fft, wu) 


being a a sufficiently large positive parameter. 

For the subsequent analysis, it is instrumental to introduce the following auxiliary 
problem: find cp e V such that 

(3.4) b{<f,w) = {u-P’^ u,w)hi(q) 
for all w G where 6(-, •) is the bilinear form 

(3.5) b{v^ w)=j‘^a^{v^ w) + r{u] w) + a{v^ w). 


We assume the validity of the following regularity result (see, e.g., [22j Theorem A.l] 
for a proof in the case a rectangular domain O). 

Assumption 3.2. Let (p be the solution of ( |3.4| ). Then it holds 

(3.6) ||^||/f3(o) < Ccl\\u — P^u\\h^(^q^) 


with P^u be the elliptie projeetion defined in (3.2) and where Cq is a positive eonstant 
only depending on Li. 

We now collect some technical results that will be useful to prove the main result 
(Theorem 3.6). 


Lemma 3.2. Let u be the solution to (2.3) and P^u be the elliptie projeetion 
defined in (3.2). Then it holds 


(3.7) ||M-P'‘M||H2(n) < 

(3.8) \\u — P^u\\hi{q) ^ 


Proof. It is worth observing that the solution u to ( |2.3| ) satifies 
(3.9) b{u, tph) = - V ■ {4>'{u)\/u) + au, tph) 

for all iph^W^. 


We first prove (3.7). Let uj G W? be a generic element to be made precise later. 


We preliminary remark that, using P^u — ui E together with Lemma 
choosing a sufficiently large, we obtain 


2.3 


and 


(3.10) bh{P^u - ui,P^u - ui) > \\P^u - 
Moreover, employing ( |3.2[ ) and ( |3.9[ ) yields 

(3.11) bhiP^u^-iPh) = {F,i’h) = b{u,i^h) G 









11 


with — V • {(j)'{u)Vu) + au. 

Thus, using ( 3.11| ) and letting be a discontinuous piecewise quadratic polyno¬ 
mial, we get 


bh{P^u — P^u — uj) = bh{P^u^ P^u — uj) — bh{ui^ P^u — uj) 

= b{u^ P^u — uj) — bh{uT^^ P^u — uj) + bh{uT^ — uj^ P^u — uj) 
= b{u, P^u — uj) — b{uT^^ P^u — uj) + bh{uj^ — uj^ P^u — uj) 


where in the last equality we apply the consistency result contained in Lemma to 
the bilinear form 


b{v,w)= y] -f^aE{v,w) + a^{v,w) + aa%{v,w). 


Eedh 


From the above identity, using ( |3.1Q ) we get 


\\P^U Ui\\‘h2(^q^ < 


(3.12) 


b{u^ P^u — ui) — b{u^ P^u — ui) + b{u — P^u — uj) 
+ hh{u^ - ui, P^u - ui). 


Let us now estimate each term on the right hand side of (3.12). From the definitions 
of the bilinear forms 6(., •) and 5(', employing the interpolation estimates given 

in Proposition , we obtain 

b{u,P^u-ui)-b{u,P^u-ui) = / ( 9 ^'(m)-7M|e)Vu-V(P'* u-w/)(ia; 

Ee^h 

(3.13) %h\\P^U - Uj\\H'i(Q,y 
Moreover, choosing uj and such that (see Proposition E3 

(3.14) \\u - ujWh^^e) + 11'^ - UttWh^^e) ^ h 

and employing the continuity properties of 5(., •) and bh{'r) we get 

(3.15) b{u - u^,P^u - ui) + bh{u^ - uj.P^u - uj) < h\\P^u - ui\\h‘^^q). 

Substituting (3.13) and ( 3.15[ ) in ( |3.12 ) and using triangle inequality together with 
(|3.14[) we get (|3.7[). 


We now prove (3.8). Taking w = u — P^u in ( |3.4[ ) yields 

(3.16) \\u — = b{ip, u — P^u) = b{ip — ipi, u — P^u) + b{(pi, u — P^u) 

We now estimate each term on the right hand side of the above equation. Choosing, 
accordingly to Proposition (3.1), (fi such that \\^ — ^ using (3.6) and 

employing the continuity property of the bilinear form 6(., •) together with (3.7) we 
get 


(3.17) b{(p - (pi,u- P’^u) < \\(p - y:’i\\H 2 (n)\\u - P'^u\\H 2 (n) < h‘^\\u - P'^u\\Hi{n). 

Using ( |3.9| ) with ipi € Wj^ we get b{ipi,u) = bh{y^i, P^u) which implies 

b{<fii, u - P^u) = bhi<fii, P^u) - b{ipi, P^u) 

= - a^{(pi,P^u)) +fh{u; (pi,P^u) - r{u; (pi,P^u) 

= ■J^Ai+A2. 
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Let and be piecewise discontinuous quadratic polynomials such that \\u — 
'^ 7 r\\H^{E) ^ hand \\(p—(Pn\\H^ (E) ^ h. Applying twice the consistency result contained 
in Lemma [2^ together with (3.7) we obtain 


EeQh EeQh 

^ (11^ - ^i\\h^(E) + 11^^ - ^7t\\h^(E)){\\P^U - u\\h^(E) + 11'^ - '^7r||i^2(£;)) 

(3-18) < h‘^\p^\H^{Q)\u\H3{Q) ^ h^||rt — P^u\\h^{q)‘ 


Let us now estimate the term A 2 . Using the definitions of r('; •, •) and •, •) we get 

^2 = W{E\e Ke{^i, P^u) - P^u)) + f (fE)\E - -/>'(«)) ■ VP^U 

Eefih 

= ■ ^ 2,1 + ^ 2 , 2 - 


Proceeding as in the bound of Ai 
if u we obtain 


and employing assumption (2.4) on the regularity 


(3.19) 


^ 2,1 ^ h^||rt - P^u\\H^(n)- 


Finally, we estimate the term ^ 2 , 2 - By employing the orthogonality property of 
projectors and denoting by (•) the projection of (•) on constants we get 


^ 2 , 2 = Y1 / (</>'(w)|£;-</>'(«)) (V<^ 7 -VP'*W-V(^-Vq dx 

Ee^h 

= E / -</>'(«)) da: 

-Eefifc 


(3.20) 


+ 


J (</>'(w)|£;-<(>'(«)) V<^-(VP'*M-Vu) da;. 


Using the interpolation estimates given in Proposition (3.1), and employing (3.6) and 
(|3.7|) together with the following inequalities 


\\\/ipi - ^(p\\l‘^(e) < ||V(^/ - ^(p\\l‘^(e) + ||V(/? - ^^\\l^(e) ^ h\\ip\\H2(^E) 

\\S/P^u\\l2^e) < llVP^rt - Vu\\l2^e) + \\"^'^\\l‘^{e) ^ (1 + h)\\u\\H2(E) 

||V(/:?||l2(£;) < ||V(/?7 - VcpWl^^e) + \\"^^\\L^iE) ^ (1 + h)||(/:?|| 772 (£;) 

||VP'*M - Vu\\l-2(e) = 1|V(P'‘m -u) + (Vm - Vm)||l2(£;) < ||P'*M - u\\h2(e) + h\\u\\m(E) 

we obtain 


(3.21) 


^ 2,2 < h?\\P'^u - u\\hi(q,). 


Combining (3.18) 
Lemma 


( |3.19D , ( |3.2lD , ( |3.17D with (|3.16D we obtain (|^. □ 

(2.3) and P^u be the elliptie projeetion 


3.3. Let u he the solution to 
defined in (3.2). Then it holds 


\\ut - (P'‘M)(||//2(n) < h 

\\ut-{P^u)t\\H^(a)<h'^- 


(3.22) 

(3.23) 
























13 


Proof. It is sufficient to observe that it holds 

bh{{P'^u)t, ^h) = Kuu ^h) + {(!>'' {u)UtWu, ^ dtWiu))\EolAP’'^^ ^h) 

Ee^h 


for all fjh G W^. Then proceeding as in Lemma 3.2 and using 

\\9t{(t>'{u))\E - 4>"{u)ut\\Loo(E) = \\6u^ - 6uUt\\Loo(E) 

we obtain the thesis. □ 


Lemma 3.4. Let u he the solution to (2.3) and P^u he the elliptie projeetion 


defined in ( |3.2[ ). Then, setting p = u — P^u and 0 = P^u — Uh, it holds 
(3.24) rh{uh,Uh;0)-fh{u;P^u,0) < \0 \h^q) (||^||l 2 (o) + ||p||l 2 (o) + + h‘^) • 


Proof. We preliminary observe that using Lemma [3^ and [3^ and proceeding as 
in the proof of [22l (3.2c)] yield P^u G lT^’^(f2), with norm bounded uniformly in 
time. Moreover, it holds 

rh{uh\Uh, 0) - fh{u\ P^u, 0) = rh{uh\P'^u, 0) - fh{u] P^u, 0) + rh{uh\Uh - P^u, 0) 

= T.{<t>'(.Uh)-(t)'{u))\EaJ^E{Eu,0) 

E^Vlh 

+ rh{uh\Uh - P’^u,0) 

= A + B. 

Let us first estimate the term A which can be written as follows 
A= Yi i^fTuh) - Apu) + fpu) - 4>'{P>^u) + 4>'{P>^u) - 

E^Vlh 

= Y {I + II + III)\EalE{P^u,0). 

E^Qj}i 

Using Lemma we obtain 

^ < E (1^1 + + \in)\E\\P'^n\\H^^E)\e\m(E) 


Ee^h 


< 


||P^rt||wuoo(Q) E \E\^^MI\ + \II\ + \III\)\e\0\hhe) 


E^Vth 


< ||P'*«||M.i,~(n) ( Y + II" + 

XEenh / 

(3-25) < \\P^u\\wr^{n){Ai + An + ^///)|6>|i^i(o). 

/ \l/2 — 

where ^(.)= (I^KOfs) • Using the definition of (•) and Lemma 


obtain 


2.2 


we 


I — ^(a^ rt/i) — £;(P^rt, P^rt)) 


= J^i(^h,E{Uh - P^U,Uh P P^U)) 


(3.26) 


^ “ ^^'^IIl2(e)(II'^/i||l2(e) + II^^'^IIl2(e)) 
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which implies 


2 \ 


E^Qh 


1/2 


(3.27) 


< 


< 




\ 1/2 


KE^Qh 


II^IU 2 (n) 


where in the last step we employed Assumption |3.1| on the regularity of Uh • 
Similarly, using the definition of (•) we have 




which yields 

Am < (||■w;l||z,~(n) + ||P'*'u||Loo(n)) ( \\P^u - M||i 2 (E)) < \\P^u - u\\l 2 ^q). 


\E^Qh 


Finally, employing q G IP 2 (E) together with Lemma 2.2 and the interpolation esti¬ 
mates of Proposition |3.1[ it is easy to prove that the following holds 

p {alE{P\ P’^u) - {P\ P>^u)) = ^ {a^EiP^u - q, P^u) - {P’^u - q, P'^u)) 


< 


< 


|^||P'*M - q\\L^{E)\\P^u\\L^(E) % |-^ {\\P^U - u\\e'2(E) + \\u “ q\\L^(E)) \\P^u\\l2(e) 


E 


— {\\P^U - u\\l^(e) + h"^) \\P'^u\\l'^(E) 


\E. 

which implies 


1/2 


(3.28) 


Aii<\\P^uh^^a) ( ( E 

EeQh 

< \\P’^u-u\\L2(Q) + h^. 


lL 2 (E) 


Employing the above estimates for A/, A// and Am into (3.25) and recalling that 
||P^i/||ivi’Oo(o) is uniformly bounded in time, we get 


(3.29) 


A < (||^||L=(n) + WP’^u - w|U2(n) + h^) \e\H^ 


(O). 


To conclude it is sufficient to estimate B. Using the definition of rh{'rr) together 


with Lemma 2.3 and Assumption |3.1| we have 


(3.30) 


B < 








□ 

Lemma 3.5. Let Vh G lU^ and e > 0. Then there exists a eonstant depending 
on e sueh that it holds 


(3.31) 
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Proof. It is straightforward to observe that it holds 


\vh\ 


H^(Q) 


AvhVh dx 


(3.32) 


= J2\^~J^^'^hVhdx^ < e\\Avh\\l2^a) + Ce\\vh\\l2^a) 


} 


where we used Cauchy-Schwarz inequality and the fact that C □ 

We are now ready to prove the following convergence result. 

Theorem 3.6. Let u he the solution to (2.3) and Uh the solution to ( |2J4l ). Then 
for all t G [0, T] it holds 


(3.33) 


Wu-Ukh^Q) 


Proof. As usual, the argument is based on the following error decomposition 
(3.34) u — Uh = {u — P^u) + {P^u — Uh) ='. p P 0. 


3.2 


In view of Lemma 
first observe that it holds 


we only need to estimate \\0\\L2(^Qy Proceeding as in [22], we 


alidt, Xh) + Xh) = aliiP'^u - Uh)t, Xh) + pa^{P’^u - uh, Xh) 

= Xh) + pa^{P’^u, Xh) 

-[ah{iuh)t, Xh) + pa^{uh, Xh)] 

= al{{P^u)t, Xh) + pa^iP^u, Xh) + rniuh, Uh; Xh)- 

Using ( |3.3[ ) and ( |3.2[ ) it holds 

pa^iP'^u, Xh) = bhiP'^u, Xh) - fh(u; P’^u, Xh) - a{P'^u, Xh) 

= {pA^u - V • {(l>'{u)Wu) + au, Xh) - rhiu; P^u, Xh) - a{P^u, Xh) 
= {pA'^u - V ■ {4)'{u)Vu), Xh) - rh{u; P'^u, Xh) + a{p, Xh)- 

Thus, we have 

alidt, Xh) + PPid, Xh) = al{{P^u)t, Xh) + {pA^u - V • {(p'{u)Vu),Xh) 

+ rh{uh, Uh\Xh) - rh{u-, P^u, Xh) + oi{p, Xh) 

= -aliPt, Xh) + {ut + x^A'^u - V ■ {4)'{u)Vu), Xh) 

+ rh{uh, Uh', Xh) - rh(u; P'^u, Xh) + ot{p, Xh) 

= a{p, Xh) - a-hipt, Xh) + rhiuh, Uh; Xh) - rhiu; P^u, Xh)- 


Taking in the above equality we get 

i3.35)alidt,0) + pa^i0,d) = aip,0) - alipt,0) + rh{uh,Uh;0) - fhiu;P’^u,9) 


which, combined with the stability properties of a^(-, •) and n^(', •) (see Lemma 2.3), 
implies the following crucial inequality 

^ (aWphHQ) + \\Pt\\LHQ))\\0\\LHQ) +rhiuh,Uh;0) -fh{u;P'^u,e). 
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Employing Lemmas |3.2| |3.3[ |3.4| and |3.5| we obtain 


(3.36) 


1 d 




which, combined with Gronwall’s lemma, yields the required estimate for ||^||l 2 (o)- 

□ 


4. Numerical results. The time discretization is performed by the Backward 
Euler method. The resulting non-linear system (2.15) at each time step is solved by 
the Newton method, using the norm of the relative residual as a stopping criterion. 
The tolerance for convergence is le — 6 . Eor the simulations, we have used a Matlab 
code and a EortranQO parallel code based on the PETSc library. The parallel tests 
were run on the EERMI linux cluster of the CINECA consortium (www.cineca.it). 


Table 1 

Test 1: and L‘^ errors and convergenee rates a eomputed on four quadrilateral meshes 

diseretizing the unit square. 


h 


a 

~ u\h^(Q) 

a 

\\uh - u\\L2(a) 

a 

1/16 

1.35e-l 

- 

8.57e-2 

- 

8.65e-2 

- 

1/32 

5.86e-2 

1.20 

2 . 20 e -2 

1.96 

2 . 20 e -2 

1.97 

1/64 

2.79e-2 

1.07 

5.53e-3 

1.99 

5.52e-3 

1.99 

1/128 

1.38e-2 

1.02 

1.37e-3 

2.01 

1.37e-3 

2.01 


4.1. Test 1: convergence to exact solution. In this test, we study the con¬ 
vergence of our VEM discretization applied to the Cahn-Hilliard equation with a 
forcing term / obtained imposing as exact solution u{x,y^t) = t cos(27rT) cos(27r^). 
The parameter 7 is set to 1/10 and the time step size At is le — 7. The and 

^2 computed at t = 0.1 on four quadrilateral meshes discretizing the unit 

square. 

The results reported in Table [T] show that in th e norm the VEM method 


In the and 


the 


converges with order 2, as predicted by Theorem 3.6 
method converges with order 1 and 2 respectively, as can be expected according to 
the EEM theory and the approximation properties of the adopted virtual space. 


4.2. Test 2: evolution of an ellipse. In this test, we consider the Cahn- 
Hilliard equation on the unit square with 7 = 1/100. The time step size At is 5e — 5. 
The initial datum uq is a piecewise constant function whose jump-set is an ellipse: 


uo{x,y) 


0.95 if 9(:c-0.5)2+ (y-0.5)2 < 1/9, 
—0.95 otherwise . 


Both a structured quadrilateral mesh and an unstructured triangular mesh (generated 
with the mesh generator of the Matlab PDEToolbox) are considered, with 49923 and 
13167 dof, respectively. As expected the initial datum uq with the ellipse-shaped 
jump-set evolves to a steady state exhibiting a circular interface; see Eigs |l(a)| and 
|l(b)| Thereafter, no motion will occur as the interface has constant curvature. 


4.3. Test 3: evolution of a cross. We use here the same domain and the 
parameters as in Test 2. The initial datum uq is a piecewise constant function whose 
jump-set has the shape of a cross; see Eigs. 3(a)[ |3(b)| and |3(c)| (t = 0). The same 
quadrilateral and triangular meshes of the Test 2 are considered, with 49923 and 13167 
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(a) Quadrilateral mesh of 16384 = 128 X 128 elements (49923 degrees of freedom). 



(b) Triangular mesh of 8576 elements (13167 degrees of freedom). 

Figure 1. Test 2: evolution of an ellipse at three temporal frames (t = 0,0.5, 1 ). 



Figure 2. Examples of Voronoi polygonal meshes (quadrilaterals, pentagons, hexagons) with 
10 (left) and 100 (right) elements. 


dof, respectively, and a Voronoi polygonal mesh (including quadrilaterals, pentagons 
and hexagons, see Fig. [^as example) with 59490 dof. As in the ellipse example, the 
initial datum uq with a cross-shaped jump-set evolves to a steady state exhibiting a 
circular interface, see Figs. |3(a')| |3(b)| and 3(c)| 


4.4. Test 4: spinoidal decomposition. Spinodal decomposition is a physical 
phenomenon consisting of the separation of a mixture of two or more components to 
bulk regions of each. It occurs when a high-temperature mixture of different compo¬ 
nents is rapidly cooled. To model this separation the initial datum uq is chosen to be 
a uniformly distributed random perturbation between -1 and 1, see Figs. 4(a)[ |4(1^ 
4(c)| (t=0). The same parameters as in Test 2 are used. We remark that the three 
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(a) Quadrilateral mesh of 16384 = 128 X 128 elements (49923 degrees of freedom). 



(b) Triangular mesh of 8576 elements discretizing the unit square (13167 degrees of freedom). 



(c) Voronoi polygonal mesh (quadrilaterals, pentagons, hexagons) of 10000 elements (59490 degrees 
of freedom). 

Figure 3. Test 3: evolution of a eross at three temporal frames (t = 0,0.05, 1). 


initial random configurations are different. We consider a quadrilateral mesh with 
49923 dof (Fig. |4(a)] ), a triangular mesh with 13167 dof (Fig. |4(b)| and a polygonal 
mesh with 59590 dof (Fig. |4(c)[ ). The separation of the two components into bulk 
regions can be appreciated quite early, see Figs. 4(a)[|4(b')||4(c)| (t=Q.Ql). This initial 
separation happens over a very small time-scale compared to the motion thereafter. 
Then, the bulk regions begin to move more slowly, and separation will continue until 
the interfaces develop a constant curvature. In the quadrilateral (Fig. |4(a)] ) and trian¬ 
gular (Fig. |4(b')] ) mesh cases, the final equilibrium configuration is the square divided 
into two rectangles, while in the polygonal (Fig. |4(c)[ ) mesh case the final equilib¬ 
rium configuration is clearly a circle. The fact that different final configurations are 
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(b) Triangular mesh of 8576 elements (13167 degrees of freedom). 



(c) Voronoi polygonal mesh (quadrilaterals, pentagons, hexagons) of 10000 elements (59490 degrees 
of freedom). 

Figure 4. Test 4- spinoidal decomposition at three temporal frames (t = 0.01,0.05,5 for the 
quadrilateral and Voronoi polygonal meshes^ t = 0.075,0.25,1.25 for the triangular mesh). 


obtained starting from different initial random configurations is consistent with the 
results in [29] . 
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